Progetto Esame DATA SCIENCE

1 Introduzione

Il Caenorhabditis elegans (C. elegans) non è un semplice verme: è una pietra miliare della neuroscienza computazionale. È l’unico organismo di cui abbiamo mappato l’intero sistema nervoso, noto come Connettoma.

Questo progetto si propone di analizzare la rete neurale del C. elegans non come un semplice oggetto biologico, ma come un Grafo Diretto e Pesato, utilizzando gli strumenti della Network Science. Il dataset utilizzato, proveniente dal progetto OpenWorm, ci permette di distinguere tra neuroni sensoriali (Input), interneuroni (Elaborazione) e neuroni motori (Output), offrendo un’opportunità unica per studiare come una struttura topologica complessa possa generare comportamenti biologici funzionali.

I dati derivano dall’aggiornamento del dataset degli anni ’80, completandone la struttura.

2 Descrizione dei Dataset e Struttura della Rete

Per questo progetto analizzeremo il sistema nervoso del C. elegans utilizzando tre file distinti provenienti dalla repository OpenWorm. Questi dati ci permettono di ricostruire un Grafo Diretto e Pesato (\(G = (V, E, W)\)) e di classificare i nodi in base alla loro funzione biologica.

2.1 1. Il Dataset Topologico: Connectome.csv

Questo è il file principale (Core Dataset) e rappresenta la Lista degli Archi (Edge List) della rete neurale.

• Contenuto: Ogni riga definisce una connessione sinaptica tra due neuroni.
• Variabili Chiave:
  • Neuron (Source): Il neurone di partenza (Presinaptico).
  • Target (Target): Il neurone di arrivo (Postsinaptico).
  • Number of Connections (Weight): Il numero di sinapsi fisiche rilevate tra i due neuroni.
  • Neurotransmitter (Attribute): Il tipo di neurotrasmettitore (es. exc per eccitatorio, inh per inibitorio).

2.2 2. Il Dataset di Input: Sensory.csv

Questo dataset funge da Tabella degli Attributi dei Nodi per l’input. Identifica quali neuroni agiscono come recettori dell’ambiente esterno.

• Contenuto: Associa specifici neuroni a funzioni sensoriali.
• Variabili Chiave:
  • Function: Il tipo di stimolo rilevato.
  • Neuron: L’ID del neurone sensoriale.

Utilizzo: Useremo questo file per aggiungere un attributo type = "Sensory" ai nodi corrispondenti nel grafo principale.

2.3 3. Il Dataset di Output: Neurons_to_Muscles.csv

Questo dataset descrive l’interfaccia tra il sistema nervoso e il sistema muscolare.

• Contenuto: Mappa le connessioni dai neuroni ai muscoli fisici.
• Variabili Chiave:
  • Origin: Il neurone motore.
  • Muscle: Il muscolo bersaglio (nota: questo non è un neurone, ma un organo effettore).
• Interpretazione: I neuroni elencati nella colonna Origin sono i Motoneuroni. Rappresentano l’ultimo stadio dell’elaborazione neurale prima dell’azione fisica.

Utilizzo: I neuroni presenti in questo file verranno etichettati con l’attributo type = "Motor".

3 Caricamento Librerie

4 Caricamento Dataset

5 Visualizzazione Grafo Originale

edges_synapse <- connectome %>%
  select(from = Neuron, to = Target, weight = `Number of Connections`) %>%
  mutate(type = "synapse")

edges_muscle <- muscles %>%
  select(from = Origin, to = Muscle, weight = `Number of Connections`) %>%
  mutate(type = "neuromuscular")

all_edges <- bind_rows(edges_synapse, edges_muscle)

g_tidy <- as_tbl_graph(all_edges, directed = TRUE)

list_sensory <- unique(sensory$Neuron)
list_motor   <- unique(muscles$Origin)
list_muscle  <- unique(muscles$Muscle)
list_hybrid  <- intersect(list_sensory, list_motor)

g_final <- g_tidy %>%
  activate(nodes) %>%
  # Uniamo le coordinate (ora chiamate x_phys, ecc.)
  left_join(distances, by = "name") %>%
  # Assegnazione Ruoli
  mutate(role = case_when(
    name %in% list_muscle ~ "Muscle",
    name %in% list_hybrid ~ "Sensorimotor",
    name %in% list_motor ~ "Motor",
    name %in% list_sensory ~ "Sensory",
    TRUE ~ "Interneuron"
  )) 

ggraph(g_final, layout = "fr") +
  geom_edge_link(aes(color = type), alpha = 0.4) +
  geom_node_point(aes(color = role), size = 3) + 
  scale_color_manual(values = c("Sensory" = "forestgreen", 
                                "Motor" = "orange", 
                                "Sensorimotor" = "purple", 
                                "Muscle" = "red", 
                                "Interneuron" = "grey")) +
  scale_edge_color_manual(values = c("synapse" = "grey", 
                                     "neuromuscular" = "red")) +
  theme_graph() +
  labs(title = "C. elegans Connectome",
       subtitle = "Rete completa con distinzione dei neuroni Ibridi")

Nella definizione dei ruoli c’è una categoria speciale chiamata “Sensorimotor”. Questi neuroni sono neuroni ibridi definiti matematicamente come l’intersezione tra l’insieme dei neuroni sensoriali e l’insieme dei neuroni motori. In termini biologici, è un neurone che ha la capacità sia di ricevere stimoli dall’ambiente (Input) sia di comandare direttamente un muscolo (Output), bypassando la rete di elaborazione intermedia.

# Calcoliamo le componenti connesse
comps <- components(g_final, mode = "weak")

6 Spiegazione topologia della Rete

L’analisi preliminare della topologia ha evidenziato che la rete non è interamente connessa (connected graph), ma presenta una componente principale e alcuni nodi isolati. Questi nodi isolati fanno parte del sistema nervoso faringeo. Questi 20 neuroni creano un “cervello autonomo” all’interno del verme, facendo si che quest’ultimo abbia un “piccolo computer” solo per gestire la deglutazione.

list_pharyngeal <- c("I1L", "I1R", "I2L", "I2R", "I3", "I4", "I5", "I6", 
                     "M2L", "M2R", "M3L", "M3R", "M4", "M5", "MCL", "MCR", 
                     "NSML", "NSMR", "M1", "MI")

g_final_ALL <- g_tidy %>%
  activate(nodes) %>%
  left_join(distances, by = "name") %>%
  mutate(role = case_when(
    name %in% list_pharyngeal ~ "Pharyngeal", # Nuova categoria!
    name %in% list_muscle ~ "Muscle",
    name %in% list_hybrid ~ "Sensorimotor",
    name %in% list_motor ~ "Motor",
    name %in% list_sensory ~ "Sensory",
    TRUE ~ "Interneuron"
  ))

ggraph(g_final_ALL, layout = "fr") +
  
  geom_edge_link(aes(color = type), alpha = 0.3, width = 0.5) +
  
  geom_node_point(aes(fill = role), shape = 21, color = "black", stroke = 0.3, size = 4) + 
  
  scale_fill_manual(values = c("Sensory" = "forestgreen", 
                                "Motor" = "orange", 
                                "Sensorimotor" = "purple", 
                                "Muscle" = "red",
                                "Pharyngeal" = "blue",
                                "Interneuron" = "grey85")) + # Grigio più chiaro per contrasto col bordo nero
  
  scale_edge_color_manual(values = c("synapse" = "grey70", 
                                     "neuromuscular" = "#FF5555")) + # Rosso leggermente più morbido
  
  theme_graph(background = "white") + # Assicura sfondo bianco pulito
  labs(title = "C. elegans Connectome",
       subtitle = "Rete completa con distinzione dei neuroni Ibridi e sistema nervoso faringeo",
       fill = "Ruolo",  
       edge_color = "Tipo Connessione")

In questa maniera abbiamo scoperto due “computer” diversi. Il Sistema Somatico (la componente gigante) deve coordinare il movimento di tutto il corpo (complesso, gerarchico). Il Sistema Faringeo (la componente piccola) deve solo gestire il ritmo della deglutizione.

7 Analisi Comparativa

La prima domanda che mi sono posto analizzando questa tipologia di rete, una volta trovate le due componenti connesse è stata in cosa differiscono i due sistemi nervosi (Somatico e faringeo).

g_marked <- g_final_ALL %>%
  activate(nodes) %>%
  mutate(comp_id = group_components(type = "weak")) # assegno un numero uguale se esiste un percorso che collega i nodi

# Identificazione degli ID 
comp_counts <- g_marked %>%
  as_tibble() %>% # Estrae la tabella dei nodi
  count(comp_id, sort = TRUE)

id_giant <- comp_counts$comp_id[1] # Somatico
id_small <- comp_counts$comp_id[2] # Faringeo

# Grafo Somatico (Giant Component)
g_somatic <- g_marked %>%
  filter(comp_id == id_giant) %>%
  mutate(subsystem = "Somatic")

# Grafo Faringeo (Pharyngeal System)
g_pharyngeal <- g_marked %>%
  filter(comp_id == id_small) %>%
mutate(subsystem = "Pharyngeal")

# Plot Somatico
p1 <- ggraph(g_somatic, layout = "fr") + 
  geom_edge_link(aes(color = type), alpha = 0.3, width = 0.5) +
  
  geom_node_point(aes(fill = role), shape = 21, color = "black", stroke = 0.3, size = 4) + 
  
  scale_fill_manual(values = c("Sensory" = "forestgreen", 
                                "Motor" = "orange", 
                                "Sensorimotor" = "purple", 
                                "Muscle" = "red",
                                "Pharyngeal" = "blue",
                                "Interneuron" = "grey85")) + # Grigio più chiaro per contrasto col bordo nero
  
  scale_edge_color_manual(values = c("synapse" = "grey70", 
                                     "neuromuscular" = "#FF5555")) + # Rosso leggermente più morbido
  theme_graph() + 
  labs(title = "Sistema Somatico (Main)")

# Plot Faringeo
p2 <- ggraph(g_pharyngeal, layout = "kk") + # 'kk' funziona meglio per grafi piccoli
  geom_edge_link(aes(color = type), alpha = 0.3, width = 0.5) +
  
  geom_node_point(aes(fill = role), shape = 21, color = "black", stroke = 0.3, size = 4) + 
  
  scale_fill_manual(values = c("Sensory" = "forestgreen", 
                                "Motor" = "orange", 
                                "Sensorimotor" = "purple", 
                                "Muscle" = "red",
                                "Pharyngeal" = "blue",
                                "Interneuron" = "grey85")) + # Grigio più chiaro per contrasto col bordo nero
  
  scale_edge_color_manual(values = c("synapse" = "grey70", 
                                     "neuromuscular" = "#FF5555")) + # Rosso leggermente più morbido
  theme_graph() + 
  labs(title = "Sistema Faringeo (Isolated)")

print(p1)

print(p2)

7.1 Confronto Densità

densita_somatic <- edge_density(g_somatic)
densita_pharyngeal <- edge_density(g_pharyngeal)

La densità del sistema somatico è 0.0198

La densità del sistema faringeo è 0.2237

Il problema è che in questa maniera sto confrontando la densità grezza di una rete di 300 nodi con una di 20. Sarebbe come barare, perché la densità tende naturalmente a scendere all’aumentare della dimensione.

Cosi ho deciso di simulare un campionamento: “Se prendessi a caso un pezzetto di cervello somatico grande quanto la faringe (20 nodi), quanto sarebbe denso?”.

L’ho eseguito 1000 volte, cosi facendo ho una prova statisticamente significativa, e non è solo un effetto della dimensione.

il P-value è un P-Value Emprico basato su ricampionamento. ho calcolato la frequenza relativa con cui questi sottografi casuali raggiungevano o superavano la densità osservata nel sistema faringeo.

Z-Score ci dice quanto è distante la nostra osservazione reale dal comportamento “medio” o “casuale”. In questo caso è stato usato per standardizzare la differenza tra la densità osservata e quella simulata. Un Z-Score positivo molto alto (es. > 3) indica che il sistema faringeo si trova nella ‘coda’ estrema destra della distribuzione nulla: è molto più denso di quanto il modello casuale possa prevedere.

densita_osservata_faringe <- edge_density(g_pharyngeal)
n_nodi_faringe <- vcount(g_pharyngeal)

# Ricampionamento
set.seed(123) # Per rendere riproducibile l'esperimento
n_permutazioni <- 1000

# Funzione che estrae un sottografo casuale e calcola la densità
simula_densita <- function() {
  # 1. Pesca N nodi a caso dal cervello somatico
  nodi_random <- sample(V(g_somatic)$name, n_nodi_faringe)
  
  # 2. Crea il sottografo indotto (mantiene solo le connessioni tra loro)
  sub_g <- induced_subgraph(g_somatic, nodi_random)
  
  # 3. Calcola la densità
  return(edge_density(sub_g))
}

# Ripetiamo la simulazione 1000 volte
distribuzione_nulla <- replicate(n_permutazioni, simula_densita())


# Z-Score
media_nulla <- mean(distribuzione_nulla)
sd_nulla <- sd(distribuzione_nulla)
z_score <- (densita_osservata_faringe - media_nulla) / sd_nulla

# Calcolo P_VALUE, quante volte il caso ha battuto la Faringe?
p_value <- sum(distribuzione_nulla >= densita_osservata_faringe) / n_permutazioni


# Visualizzazione
df_plot <- data.frame(densita = distribuzione_nulla)

ggplot(df_plot, aes(x = densita)) +
  # Istogramma della distribuzione nulla (il "rumore di fondo")
  geom_histogram(fill = "lightgrey", color = "grey", bins = 30) +
  # Linea rossa: osservazione reale
  geom_vline(xintercept = densita_osservata_faringe, color = "red", linetype = "dashed", linewidth = 1.5) +
  labs(title = "Test di Significatività della Densità",
       subtitle = paste("Faringe (Rossa) vs 1000 Sottografi Somatici Casuali (Grigio)\nP-value:", p_value),
       x = "Densità",
       y = "Frequenza") +
  theme_minimal()

7.2 Confronto Average Path Length e Transitivity

calcola_metriche <- function(grafo, nome) {
  data.frame(
    Sistema = nome,
    # Transitività (Clustering/Compattezza locale)
    Transitivity = transitivity(grafo, type = "global"),
    # Average Path Length (Efficienza di trasmissione)
    AvgPathLength = mean_distance(grafo, directed = TRUE) #aggiungo directed = true per rispettare il flusso
  )
}

df_somatic <- calcola_metriche(g_somatic, "Somatico (Main)")
df_pharyngeal <- calcola_metriche(g_pharyngeal, "Faringeo (Isolated)")


df_compare <- bind_rows(df_somatic, df_pharyngeal)

L’alta transitività della faringe suggerisce una struttura a ‘Clique’, dove i neuroni tendono a formare gruppi chiusi e interconnessi, tipico dei moduli funzionali autonomi.

Per il Percorso Breve invece lo analizziamo come per la densità, non possiamo confrontare direttamente i numeri perché una rete piccola ha naturalmente percorsi più brevi. Dobbiamo chiedere: “Il sistema Faringeo è più efficiente di un pezzo casuale di cervello Somatico?”

set.seed(123)
n_permutazioni <- 1000

# Calcoliamo i valori reali
path_osservato_faringe <- mean_distance(g_pharyngeal, directed = TRUE)
n_nodi_faringe <- vcount(g_pharyngeal)

# Funzione di simulazione
simula_path <- function() {
  # 1. estrazione nodi casuali dal somatico
  nodi_random <- sample(V(g_somatic)$name, n_nodi_faringe)
  
  # 2. Crea il sottografo
  sub_g <- induced_subgraph(g_somatic, nodi_random)
  
  # 3. Gestione disconnessione
  comps <- components(sub_g)
  giant_sub <- induced_subgraph(sub_g, which(comps$membership == which.max(comps$csize)))
  
  # Se troppo piccolo, NA
  if(vcount(giant_sub) < 2) return(NA)
  
  return(mean_distance(giant_sub, directed = TRUE))
}

# Eseguiamo il bootstrap
distribuzione_path <- replicate(n_permutazioni, simula_path())
distribuzione_path <- distribuzione_path[!is.na(distribuzione_path)]

# Calcolo P-Value
p_value_path <- sum(distribuzione_path <= path_osservato_faringe) / length(distribuzione_path)

trans_osservata_faringe <- transitivity(g_pharyngeal, type = "global")
n_nodi_faringe <- vcount(g_pharyngeal)

# setto il campione per riproducibilità
set.seed(123)
n_permutazioni <- 1000

simula_transitivity <- function() {
  #Pesca nodi casuali dal somatico
  nodi_random <- sample(V(g_somatic)$name, n_nodi_faringe)
  
  #Crea sottografo
  sub_g <- induced_subgraph(g_somatic, nodi_random)
  
  #Calcola Transitività
  # grafo non ha connessioni, transitivity restituisce NaN.
  t_val <- transitivity(sub_g, type = "global")
  
  # Se NaN (nessuna triade), lo consideriamo 0
  if(is.nan(t_val)) return(0)
  
  return(t_val)
}

# Esecuzione
distribuzione_trans <- replicate(n_permutazioni, simula_transitivity())


# P-VALUE
p_value_trans <- sum(distribuzione_trans >= trans_osservata_faringe) / n_permutazioni

# Z-Score
media_nulla <- mean(distribuzione_trans)
sd_nulla <- sd(distribuzione_trans)
z_score <- (trans_osservata_faringe - media_nulla) / sd_nulla

df_plot <- data.frame(valore = distribuzione_trans)

ggplot(df_plot, aes(x = valore)) +
  geom_histogram(fill = "lightgrey", color = "grey", bins = 30) +
  geom_vline(xintercept = trans_osservata_faringe, color = "blue", linetype = "dashed", size = 1.5) +
  labs(title = "Test di Significatività: Transitività",
       subtitle = paste("Faringe (Blu) vs Random Somatico. P-val:", p_value_trans),
       x = "Transitività Globale", y = "Frequenza") +
  theme_minimal()

7.3 Conclusioni Comparazioni

dati_plot <- data.frame(
  Sistema = c(rep("Somatico (Main)", 3), rep("Faringeo (Isolated)", 3)),
  Metrica = rep(c("Density", "Transitivity", "AvgPathLength"), 2),
  Valore = c(0.02, 0.10, 3.5,   # Valori approssimativi Somatici
             0.23, 0.40, 2.8)   # Valori approssimativi Faringei (dai tuoi risultati)
)

ggplot(dati_plot, aes(x = Sistema, y = Valore, fill = Sistema)) +
  geom_bar(stat = "identity", position = "dodge", width = 0.6, color = "black") +
  facet_wrap(~Metrica, scales = "free_y") +
  scale_fill_manual(values = c("Somatico (Main)" = "#4E79A7", 
                               "Faringeo (Isolated)" = "#E15759")) +
  labs(title = "Confronto Topologico: Somatico vs Faringeo",
       subtitle = "La Faringe si distingue per Densità e Transitività (Robustezza), non per Efficienza.",
       y = "Valore Metrica", x = "") +
  theme_minimal() +
  theme(legend.position = "none",
        strip.text = element_text(face = "bold", size = 12),
        plot.title = element_text(face = "bold")) +
  geom_text(aes(label = round(Valore, 2)), vjust = -0.5)

L’analisi topologica comparativa, supportata da test di ricampionamento (Bootstrap, n=1000), ha evidenziato differenze strutturali critiche tra la Componente Gigante (Sistema Somatico) e la componente isolata (Sistema Faringeo), riflettendo le diverse specializzazioni funzionali dei due circuiti.

• Densità Reale Faringe: 0.2237
• Densità Media Somatico (Simulato): 0.02
• Z-Score: 20.08
• P-Value: 0
  • Densità (Connettività Locale): Posso concludere che nonostante la dipendenza della densità dalla dimensione della rete (\(N\)), il test di permutazione (n=1000) conferma che l’alta densità del sistema faringeo è una proprietà topologica intrinseca e statisticamente significativa (\(p < 0.001, Z-score > 3\)), e non un mero artefatto di scala. ideale per attivare massivamente e simultaneamente i muscoli della deglutizione (comportamento “tutto o niente”).
• Path Length Faringe: 2.797
• Path Length Random (Media): 3.412
• P-Value: 0.5390782
  • Average Path Length (Efficienza Globale): Sebbene il sistema faringeo sia molto efficiente (\(L \approx 2.8\)), la differenza rispetto al tessuto somatico casuale non è statisticamente significativa (P \(\approx\) 0.54). la Small-Worldness non è una prerogativa esclusiva della faringe, ma è una proprietà ubiquitaria dell’intero connettoma del C. elegans. L’evoluzione ha ottimizzato l’intero cervello per minimizzare la distanza di comunicazione, indipendentemente dal modulo funzionale considerato.
• Transitività Faringe (Reale): 0.3976
• Transitività Random (Media): 0.1003
• Z-Score: 1.49
• P-Value: 0.091
  • Transitività (Robustezza e Ridondanza):Il sistema faringeo presenta un coefficiente di clustering circa 4 volte superiore al caso (\(0.40\) vs \(0.10\)). Sebbene il p-value sia marginale (P \(\approx\) 0.09) a causa dell’alta varianza campionaria su reti piccole, l’Effect Size è marcato. il sistema faringeo tende a organizzarsi in clique e triangoli chiusi.

8 Analisi Rete Somatica

Procediamo l’analisi con la rete somatica, che ci permette di fare analisi più approfondite essendo più grande.

8.1 Analisi Distribuzioni dei gradi

# Estrazione dei gradi In/out/total della componente gigante
dati_gradi <- data.frame(
  Nome = V(g_somatic)$name,
  Total = degree(g_somatic, mode = "all"),
  In = degree(g_somatic, mode = "in"),
  Out = degree(g_somatic, mode = "out")
)

# Visualizzazione 
dati_plot <- dati_gradi %>%
  pivot_longer(cols = c(Total, In, Out), names_to = "Tipo", values_to = "Grado") %>%
  mutate(Tipo = factor(Tipo, levels = c("In", "Out", "Total"))) # Ordine logico

ggplot(dati_plot, aes(x = Grado)) +
  geom_histogram(fill = "steelblue", color = "black", binwidth = 2) +
  facet_wrap(~Tipo, scales = "free") +
  labs(title = "Distribuzione dei Gradi (Sistema Somatico)",
       subtitle = "Notare la coda lunga (Long Tail) in tutte le metriche",
       y = "Frequenza (Nodi)", x = "Grado (k)") +
  theme_minimal()

#vettore per i colori 
my_node_colors <- c("Sensory" = "forestgreen", 
                    "Motor" = "orange", 
                    "Sensorimotor" = "purple", 
                    "Muscle" = "red", 
                    "Pharyngeal" = "blue", 
                    "Interneuron" = "grey85")

my_edge_colors <- c("synapse" = "grey70", 
                    "neuromuscular" = "#FF5555")

# Calcoliamo i gradi e fissiamo il layout 
g_plot <- g_somatic %>%
  activate(nodes) %>%
  mutate(
    deg_total = centrality_degree(mode = "all"),
    deg_in    = centrality_degree(mode = "in"),
    deg_out   = centrality_degree(mode = "out")
  )

set.seed(123) # Per riproducibilità
# Calcoliamo il layout fisso (es. Fruchterman-Reingold)
layout_fixed <- create_layout(g_plot, layout = "fr") 


plot_colored_network <- function(layout_data, size_metric, title_text) {
  
  ggraph(layout_data) +
  
    geom_edge_link(aes(color = type), 
                   alpha = 0.4, 
                   width = 0.6,
                   arrow = arrow(length = unit(2, 'mm'), type = "closed"),
                   # end_cap assicura che la freccia non finisca SOTTO il bordo nero del nodo
                   end_cap = circle(4, 'mm')) + 
    
    geom_node_point(aes(size = {{size_metric}}, fill = role), 
                    shape = 21, 
                    color = "black", 
                    stroke = 0.4, 
                    alpha = 0.9) +
    
    scale_fill_manual(values = my_node_colors, name = "Ruolo Neurale") +
    
    scale_edge_color_manual(values = my_edge_colors, name = "Tipo Connessione") +
    
    scale_size_continuous(range = c(2, 10), limits = c(0, NA), name = "Grado") +
    
    theme_graph(base_family = "sans", background = "white") + # Forza sfondo bianco pulito
    labs(title = title_text) +
    theme(
      plot.title = element_text(face = "bold", hjust = 0.5, size = 16),
      legend.position = "right" # Mantiene la legenda ordinata a destra
    )
}

# CREAZIONE DEI TRE GRAFICI

p1 <- plot_colored_network(layout_fixed, deg_total, "1. Total Degree\n(Importanza Globale)")
p2 <- plot_colored_network(layout_fixed, deg_in,    "2. In-Degree\n(Input / Integrazione)")
p3 <- plot_colored_network(layout_fixed, deg_out,   "3. Out-Degree\n(Output / Broadcasting)")

# COMPOSIZIONE FINALE (PATCHWORK)

final_plot <- (p1 + p2 + p3) + 
  plot_layout(guides = "collect") + # Raccoglie la legenda a destra
  plot_annotation(
    title = 'Analisi grado del Sistema Somatico',
    subtitle = 'Confronto della centralità (dimensione nodi) mantenendo i ruoli funzionali (colore)',
    theme = theme(plot.title = element_text(size = 16, face = "bold"))
  ) & theme(legend.position = "bottom") # Sposta la legenda in basso

# Visualizza
print(final_plot)

L’analisi visiva delle tre metriche di centralità

• Out-Degree (Broadcasting / Divergenza):
  • Significato Questi nodi agiscono come Broadcasters. Il loro compito è trasdurre uno stimolo ambientale (es. un tocco sul naso) e disseminare l’informazione (“urlarla”) verso il maggior numero possibile di nodi.
• In-Degree (Integrazione / Convergenza):
  • Significato Questi nodi fungono da Integratori. Per eseguire un’azione critica come la contrazione muscolare, non possono affidarsi a un singolo input; necessitano di ricevere “l’ordine finale” convergente da molteplici nodi. Questa ridondanza in entrata serve a filtrare il rumore e prevenire attivazioni motorie accidentali.
• Total Degree (Hub Centrali / Routing):
  • Nodi Dominanti: I nodi con il grado totale più elevato sono tipicamente gli Interneuroni di Comando.
  • Significato Questi nodi rappresentano il nucleo computazionale (“Core”) della rete. Posizionati strategicamente tra sensori e motori, possiedono sia un alto In-Degree (raccolta dati) che un alto Out-Degree (invio comandi), gestendo la logica complessa del comportamento (es. inversione di marcia).

stats_by_role <- g_somatic %>%
  activate(nodes) %>%
  mutate(
    d_out = centrality_degree(mode = "out"),
    d_in = centrality_degree(mode = "in")
  ) %>%
  as_tibble() %>%
  # media
  group_by(role) %>%
  summarise(
    N_Neuroni = n(),
    Avg_Out_Degree = mean(d_out),
    Avg_In_Degree = mean(d_in)
  ) %>%
  arrange(desc(Avg_Out_Degree))

# Stampiamo la tabella
kable(stats_by_role, digits = 2, caption = "Grado Medio per Ruolo Funzionale")

Grado Medio per Ruolo Funzionale

role	N_Neuroni	Avg_Out_Degree	Avg_In_Degree
Interneuron	88	11.38	11.74
Sensorimotor	13	10.92	4.62
Motor	105	9.09	7.74
Sensory	73	8.85	3.95
Muscle	94	0.00	5.84

L’analisi quantitativa dei gradi medi per categoria rivela il flusso dettagliato dell’informazione:

1. Origine (Sensory): I neuroni sensoriali mostrano un chiaro sbilanciamento verso l’output (\(K_{out} \approx 8.8\) vs \(K_{in} \approx 3.9\)), confermando il loro ruolo di broadcasters primari.
2. Elaborazione (Interneuron): Mantengono un bilanciamento quasi perfetto (\(K_{out} \approx 11.4\) vs \(K_{in} \approx 11.7\)), agendo come hub di ricircolo e integrazione.
3. Amplificazione (Motor): I motoneuroni mostrano un Out-Degree superiore all’In-Degree (\(9.09 > 7.74\)). Questo riflette la loro funzione di divergenza sinaptica verso i muscoli: un singolo neurone motorio deve innervare molteplici cellule muscolari per coordinare la contrazione.
4. Destinazione (Muscle): I nodi muscolari rappresentano i veri “Sinks” della rete, con un In-Degree netto (~5.8) e un Out-Degree nullo (0.00), segnando la fine del flusso informativo topologico.

8.2 Analisi Power Law e esponenziale della distribuzione dei gradi

# Dati
x <- dati_gradi$Total
x <- x[x > 0] 

# Fit Power Law eseguo la power law su un Xmin per tenere in considerazione i nodi importanti
m_pl <- displ$new(x)
est_pl <- estimate_xmin(m_pl)
m_pl$setXmin(est_pl)

# Fit Esponenziale (sullo stesso xmin)
m_exp <- disexp$new(x)
m_exp$setXmin(est_pl$xmin)
est_exp <- estimate_pars(m_exp)
m_exp$setPars(est_exp)

# Confronto 
comp <- compare_distributions(m_pl, m_exp)
test_stat <- comp$test_statistic
p_val <- comp$p_value

# Se p_val è NULL o NA, lo impostiamo ad un valore fittizio
if(is.null(p_val) || is.na(p_val)) {
  p_val_txt <- "NA (Non calcolabile)"
  p_val_num <- 1.0 # Consideriamo non significativo per prudenza
} else {
  p_val_txt <- as.character(round(p_val, 4))
  p_val_num <- p_val
}

# Estrazione parametro lambda (Esponenziale)
lambda_exp <- m_exp$getPars()


cat("Soglia ottimale (xmin):", est_pl$xmin, "\n")

## Soglia ottimale (xmin): 21

cat("Esponente Power-Law (alpha):", round(est_pl$pars, 3), "\n")

## Esponente Power-Law (alpha): 3.812

cat("Esponente Esponenziale (lambda):", round(lambda_exp, 5), "\n")

## Esponente Esponenziale (lambda): 0.09189

# Grafico del Fit
plot(m_pl, main="Fit: Power-Law (Rosso) vs Esponenziale (Blu)", xlab="Grado", ylab="CCDF")
lines(m_pl, col="red", lwd=2)
lines(m_exp, col="blue", lwd=2, lty=2)

AIC_pl <- 2 * 1 - 2 * dist_ll(m_pl)  
AIC_exp <- 2 * 1 - 2 * dist_ll(m_exp) 

cat("AIC Power-Law:", AIC_pl, "\n")

## AIC Power-Law: 550.3318

cat("AIC Esponenziale:", AIC_exp, "\n")

## AIC Esponenziale: 557.5526

delta_AIC <- AIC_pl - AIC_exp

if (delta_AIC < -2) {
  cat("Verdetto: La Power-Law è il modello migliore (AIC più basso).\n")
} else if (delta_AIC > 2) {
  cat("Verdetto: L'Esponenziale è il modello migliore (AIC più basso).\n")
} else {
  cat("Verdetto: I due modelli sono statisticamente indistinguibili (AIC simile).\n")
}

## Verdetto: La Power-Law è il modello migliore (AIC più basso).

8.2.1 Power-Law vs Distribuzione Esponenziale

L’analisi della distribuzione dei gradi confronta due modelli teorici fondamentali per descrivere la topologia della rete.

8.2.1.1 Power-Law (Legge di Potenza)

• Il Concetto: Il meccanismo sottostante è il Preferential Attachment (“Rich get richer”). I nodi che possiedono già molte connessioni hanno una probabilità maggiore di attrarne di nuove rispetto ai nodi poco connessi.
• Tipo di Rete: Si definisce Scale-Free. In queste reti non esiste un “neurone tipico” (la media non è rappresentativa). La topologia è caratterizzata da una vasta maggioranza di nodi con pochissime connessioni e da un numero esiguo di “Hub” (Mostri) con un grado teoricamente illimitato.
• Significato Biologico (Teorico): Se la rete fosse una Power-Law pura (con un esponente tipico \(2 < \alpha < 3\)), significherebbe che il sistema non è soggetto a vincoli fisici stringenti: un singolo neurone potrebbe potenzialmente connettersi all’intero cervello senza incorrere in costi energetici o spaziali proibitivi.
• Analisi del parametro \(\alpha\) (3.81): Il valore osservato è troppo alto rispetto allo standard delle reti scale-free classiche (\(\alpha \approx 2.5\)). Un esponente di 3.81 indica una pendenza molto ripida nella coda della distribuzione: la probabilità di trovare super-hub crolla molto più velocemente di quanto previsto da una legge di potenza pura.

8.2.1.2 Distribuzione Esponenziale

• Formula: \(P(k) \sim e^{-\lambda k}\).
• Il Concetto: Il meccanismo dominante è il Costo di Crescita. Ogni nuova connessione comporta un costo metabolico e un costo volumetrico.
• Tipo di Rete: Si definisce “Democratica” o Random. Esiste una scala caratteristica (un grado medio tipico) attorno alla quale si concentrano quasi tutti i nodi. I nodi con grado altissimo sono eventi statisticamente improbabili.
• Significato Biologico (Fisico): Questo modello riflette i vincoli fisici reali. Il C. elegans è un organismo microscopico con un corpo tubolare stretto. Non vi è lo spazio fisico affinché un singolo neurone possa ospitare migliaia di filamenti (assoni/dendriti). La coda della distribuzione viene quindi “tagliata” (Truncated) perché la formazione di hub giganti è fisicamente impossibile.

Ho usato Xmin, per dire di Ignorare la massa di nodi piccoli e normali che confondono i calcoli, e partire da dove inizia il comportamento speciale degli HUB, poichè La legge di potenza è un fenomeno che emerge solo quando un nodo supera una certa soglia di importanza (\(x_{min}=21\)). Sotto quella soglia, il sistema è governato da altre regole (random o geometriche).

Conclusioni:

La rete sfrutta i vantaggi degli Hub per l’efficienza comunicativa (zona scale-free), ma subisce un taglio netto nella coda estrema imposto dai limiti biofisici dell’organismo (componente esponenziale). Alla luce dei dati (\(\alpha \approx 3.8\), \(\lambda \approx 0.09\)), possiamo concludere che la rete del C. elegans non ricade perfettamente in nessuno dei due modelli estremi, ma presenta caratteristiche ibride tipiche delle reti biologiche fisiche, Truncated Power-Law. L’esponente \(\alpha\) elevato (3.8) indica che la formazione di ‘Super-Hub’ è molto più rara rispetto a una rete Scale-Free teorica. Perché tende all’Esponenziale: Il valore di \(\lambda\) (0.09) suggerisce che i costi fisici e metabolici di cablaggio agiscono come un fattore limitante (“cutoff”), impedendo ai neuroni di accumulare un numero infinito di connessioni.

Tramite il criterio AIC la power law è migliore.

Ci sono HUB per garantire l’efficienza, ma non possono accumulare connessioni all’infinito per via dell’organismo

8.3 Assortività

assortativity_val <- assortativity_degree(g_somatic, directed = TRUE)

N_TOP <- 15

# Calcoliamo i gradi
g_stats <- g_somatic %>%
  activate(nodes) %>%
  mutate(total_deg = centrality_degree(mode = "all")) %>%
  as_tibble()

# Troviamo chi sono i Top 15 (i Nomi)
top_nodes_names <- g_stats %>%
  arrange(desc(total_deg)) %>%
  slice(1:N_TOP) %>%
  pull(name)

# Creiamo un grafo fatto SOLO da questi 15 nodi (Induce Subgraph)
g_club <- induced_subgraph(as.igraph(g_somatic), vids = top_nodes_names)

# Calcoliamo la densità: 
# Quanti collegamenti ci sono tra loro RISPETTO a quelli possibili?
density_club <- edge_density(g_club)
density_global <- edge_density(g_somatic)

# Disegniamo solo questi 15 per vedere se formano una palla densa
set.seed(123)
ggraph(as_tbl_graph(g_club), layout = "fr") +
  geom_edge_link(alpha = 0.5, color = "red") +
  geom_node_point(size = 5, color = "darkred") +
  geom_node_text(aes(label = name), repel = TRUE, color = "black") +
  labs(title = paste("Il Rich Club dei Top", N_TOP),
       subtitle = paste("Densità interna:", round(density_club*100, 1), "%")) +
  theme_void()

• ANALISI RICH CLUB (Top 15 Nodi)
  • Densità Globale della Rete: 0.0198
  • Densità Interna dei Top 15: 0.3429
  • Rich Club Ratio: 17.34 x

# Calcoliamo i gradi
g_rich_ready <- g_somatic %>%
  activate(nodes) %>%
  mutate(total_deg = centrality_degree(mode = "all"))

# Troviamo il valore di taglio (il grado del 15esimo neurone)
threshold_value <- g_rich_ready %>%
  as_tibble() %>%
  arrange(desc(total_deg)) %>%
  slice(N_TOP) %>%
  pull(total_deg)

g_richviz <- g_rich_ready %>%
  activate(nodes) %>%
  mutate(
    # Creiamo un "flag" (TRUE/FALSE): Sei un VIP?
    is_rich_club = total_deg >= threshold_value,
    # Creiamo etichette SOLO per i VIP
    rich_label = ifelse(is_rich_club, name, NA),
    # Assegniamo colori specifici per il plot
    node_color = ifelse(is_rich_club, "gold", "grey80"),
    node_alpha = ifelse(is_rich_club, 1, 0.3),
    node_size  = ifelse(is_rich_club, 6, 2)
  ) %>%
  activate(edges) %>%
  mutate(
    # Un arco è "Rich" solo se connette DUE nodi Rich
    # .N() ci permette di guardare i dati dei nodi dall'interno degli archi
    from_rich = .N()$is_rich_club[from],
    to_rich = .N()$is_rich_club[to],
    is_rich_edge = from_rich & to_rich,
    # Colori e trasparenze per gli archi
    edge_color = ifelse(is_rich_edge, "red", "grey90"),
    edge_alpha = ifelse(is_rich_edge, 0.8, 0.1),
    edge_width = ifelse(is_rich_edge, 1, 0.3)
  )

set.seed(123) # Importante per mantenere la stessa forma degli altri grafici
layout_rich <- create_layout(g_richviz, layout = "fr")

ggraph(layout_rich) +
  # Disegniamo PRIMA lo sfondo (grigio chiaro)
  geom_edge_link(aes(filter = !is_rich_edge, color = edge_color, alpha = edge_alpha, width = edge_width), show.legend = FALSE) +
  # Disegniamo POI le connessioni "Rich" (rosse) sopra tutto
  geom_edge_link(aes(filter = is_rich_edge, color = edge_color, alpha = edge_alpha, width = edge_width), show.legend = FALSE) +

  geom_node_point(aes(filter = !is_rich_club, color = node_color, alpha = node_alpha, size = node_size), show.legend = FALSE) +

  geom_node_point(aes(filter = is_rich_club, color = node_color, alpha = node_alpha, size = node_size), show.legend = FALSE) +


  geom_node_text(aes(label = rich_label), repel = TRUE, size = 4, fontface = "bold", color = "black") +

  scale_edge_color_identity() + # Usa i colori esatti definiti (red/grey)
  scale_edge_alpha_identity() +
  scale_edge_width_identity() +
  scale_color_identity() +      # Usa i colori esatti definiti (gold/grey)
  scale_alpha_identity() +
  scale_size_identity() +

  theme_graph(base_family = "sans") +
  labs(title = "Visualizzazione del 'Rich Club' Core",
       subtitle = "I Top 15 Hub (Oro) e le loro interconnessioni dirette (Rosso) sono evidenziati.\nIl resto della rete è in secondo piano.",
       caption = "Nota: La densità di connessioni rosse tra i nodi dorati spiega la componente assortativa della rete.") +
  theme(plot.title = element_text(face = "bold", size = 16))

8.3.1 Analisi dell’Assortatività

Il calcolo del coefficiente di assortatività per grado restituisce un valore di \(r \approx 0.0275\). Sebbene questo valore sia numericamente prossimo allo zero (suggerendo neutralità), nel contesto specifico del connettoma di C. elegans esso riflette un bilanciamento tra due tendenze topologiche opposte:

1. Tendenza Disassortativa (Controllo): La struttura a stella tipica delle reti biologiche, dove pochi hub centrali coordinano molti nodi a basso grado, esercita una pressione verso valori negativi.
2. Tendenza Assortativa (Core Integrativo): La presenza evidenziata in precedenza di un forte “Rich Club” di interneuroni che sono densamente interconnessi tra loro, esercita una pressione verso valori positivi.

Conclusione: Il valore neutro non indica casualità, ma rivela che il sistema somatico opera con una Architettura Ibrida: possiede un nucleo centrale denso e collaborativo (Assortativo) che collettivamente gestisce una periferia gerarchica (Disassortativa).

9 Analisi Betwenness Centrality

Un ulteriore domanda a cui volevo darmi una risposta dati i risultati precedenti è la seguente: “I neuroni più connessi (Degree) sono gli stessi che controllano il flusso di informazioni (Betweenness)? O sono diversi?”

9.1 Mappa Centralità Somatica

# Calcoliamo le metriche 
g_viz_final <- g_somatic %>%
  activate(nodes) %>%
  mutate(
    # Calcoliamo le metriche
    deg_total = centrality_degree(mode = "all"),
    betw = centrality_betweenness(directed = TRUE, normalized = TRUE),
    
    # uso la radice quadrata 
    betw_viz = sqrt(betw) 
  )

# layout 
set.seed(123) 
layout_final <- create_layout(g_viz_final, layout = "fr") # Fruchterman-Reingold

ggraph(layout_final) +
 
  geom_edge_link(color = "grey85", alpha = 0.4, width = 0.2) +
  
  geom_node_point(
    aes(
      size = deg_total,   
      color = betw_viz  
    ),
    shape = 21, # Forma che supporta riempimento (color) e bordo (stroke)
    stroke = 1.5, # Spessore del bordo bianco
    fill = "white" # Colore di base se qualcosa va storto (non dovrebbe vedersi)
  ) +
  
  
  # Calibriamo per non avere nodi enormi che coprono tutto
  scale_size_continuous(
    range = c(2, 10), # Dimensione minima e massima dei pallini
    name = "Total Degree\n(Dimensione)"
  ) +
  
  #  Usiamo una palette "calda" ad alto contrasto
  scale_color_viridis_c(
    option = "plasma", 
    direction = -1, # Invertiamo: Blu=Basso, Giallo/Rosso=Alto (più intuitivo)
    name = "Betweenness\n(Colore, sqrt)"
  ) +
  
  # Top assoluti
  # Decommenta nomi dei "capi"
  # geom_node_text(aes(label = ifelse(deg_total > quantile(deg_total, 0.98), name, NA)),
                 # repel = TRUE, size = 3, fontface = "bold", color = "black") +

  theme_graph(base_family = "sans", background = "white") +
  labs(
    title = "Mappa Centralità Somatica(Degree e Betweennes)",
    subtitle = "Dimensione = Popolarità Locale (Degree) | Colore = Controllo del Flusso Globale (Betweenness)",
    caption = "Nota: I nodi grandi e con colore bluastro sono il 'Core' operativo della rete.\nLa scala colore usa la radice quadrata della betweenness per migliorare il contrasto."
  ) +
  theme(plot.title = element_text(face = "bold", size = 14))

9.2 Correlazione Degree e Betweenness

# Estraiamo i dati aggiornati dal grafo
df_centrality <- g_somatic %>%
  activate(nodes) %>%
  mutate(
    Total_Degree = centrality_degree(mode = "all"),
    # Betweenness normalizzata per confrontare reti diverse se necessario
    Betweenness  = centrality_betweenness(directed = TRUE, normalized = TRUE)
  ) %>%
  as_tibble() %>%
  select(name, role, Total_Degree, Betweenness)

# Usiamo SPEARMAN perché i dati non seguono una distribuzione normale (Power Law)
cor_val <- cor(df_centrality$Total_Degree, df_centrality$Betweenness, method = "spearman")

# Calcoliamo quanti dei Top 10 Degree sono anche nei Top 10 Betweenness
top_10_deg <- df_centrality %>% arrange(desc(Total_Degree)) %>% slice(1:10) %>% pull(name)
top_10_bet <- df_centrality %>% arrange(desc(Betweenness))  %>% slice(1:10) %>% pull(name)
overlap_count <- length(intersect(top_10_deg, top_10_bet))

ggplot(df_centrality, aes(x = Total_Degree, y = Betweenness)) +
  
  # Mappiamo il colore alla colonna 'role'
  geom_point(aes(color = role), size = 3, alpha = 0.7) +
  
  # Linea di tendenza
  geom_smooth(method = "loess", color = "grey30", linetype = "dashed", size = 0.5, se = FALSE) +
  
  # Etichette intelligenti (Top 5% o Gatekeepers se ne abbiamo trovati)
  geom_text_repel(data = filter(df_centrality, 
                                Total_Degree > quantile(Total_Degree, 0.95) | 
                                Betweenness > quantile(Betweenness, 0.95)),
                  aes(label = name), 
                  size = 3.5, box.padding = 0.4, max.overlaps = 20, 
                  color = "black", fontface = "bold") +
  
  scale_color_manual(values = my_node_colors, name = "Ruolo Funzionale") +
  
  labs(title = "Correlazione Degree vs Betweenness",
       subtitle = paste0("Spearman Rho: ", round(cor_val, 3), 
                         " | Overlap nei Top 10: ", overlap_count, "/10 neuroni"),
       x = "Total Degree (Popolarità Locale)",
       y = "Betweenness Centrality (Controllo Globale)") +
  
  theme_minimal() +
  theme(plot.title = element_text(face = "bold", size = 14),
        legend.position = "bottom",
        legend.title = element_text(face = "bold"))

## `geom_smooth()` using formula = 'y ~ x'

9.2.1 Correlazione tra Importanza Locale e Controllo Globale

Per comprendere se il controllo del flusso informativo è esercitato dagli stessi nodi che possiedono il maggior numero di connessioni, abbiamo analizzato la relazione tra Total Degree e Betweenness Centrality.

• Risultati Quantitativi:
  • Correlazione di Spearman: 0.896 (Valore molto alto indica forte sovrapposizione).
  • Overlap dei Top 10: Il 50% dei nodi presenti nella top 10 per grado sono presenti anche nella top 10 per betweenness.

Interpretazione Topologica: Il grafico di dispersione (Scatter Plot) evidenzia una forte relazione positiva. I nodi del “Rich Club” (es. AVAL, AVAR, PVCL) si posizionano nel quadrante in alto a destra (Alto Grado / Alta Betweenness). Ciò conferma che nel sistema somatico di C. elegans non vi è distinzione tra “Hub” e “Broker”: i centri di comando integrano le informazioni (funzione Hub) e agiscono contemporaneamente come autostrade obbligate per il segnale (funzione Bridge). cosi facendo la rete riesce a massimizzare la velocità di trasmissione ma concentra la vulnerabilità in pochi nodi specifici.

9.3 Analisi delle 4 centralità per Neurone

# Metriche
df_complete <- g_somatic %>%
  activate(nodes) %>%
  mutate(
    Degree      = centrality_degree(mode = "all"),
    Betweenness = centrality_betweenness(directed = TRUE, normalized = TRUE),
    Closeness   = centrality_closeness(mode = "out", normalized = TRUE), 
    PageRank    = centrality_pagerank(directed = TRUE)
  ) %>%
  as_tibble() %>%
  # Selezioniamo solo le colonne che ci servono
  select(role, Degree, Betweenness, Closeness, PageRank) %>%
  # Trasformiamo in formato lungo per il grafico
  pivot_longer(cols = -role, names_to = "Metrica", values_to = "Valore")

ggplot(df_complete, aes(x = role, y = Valore, fill = role)) +
  # Boxplot + Jitter (Punti)
  geom_boxplot(alpha = 0.7, outlier.shape = NA) +
  geom_jitter(width = 0.2, size = 0.8, alpha = 0.3) +
  
  # Dividiamo in 4 grafici indipendenti
  facet_wrap(~Metrica, scales = "free_y", nrow = 2) +
  
  # Stile e Colori
  scale_fill_manual(values = my_node_colors) +
  theme_minimal() +
  labs(title = "Profilo Funzionale Completo dei Neuroni",
       subtitle = "Confronto delle 4 centralità chiave per ruolo",
       x = NULL, y = "Valore Normalizzato") +
  theme(legend.position = "none",
        strip.text = element_text(face = "bold", size = 12),
        axis.text.x = element_text(angle = 45, hjust = 1, size = 10))

9.3.1 Analisi delle Centralità nella Rete Neurale di C. elegans

La figura mostra la distribuzione delle quattro principali misure di centralità (Betweenness, Closeness, Degree e PageRank) calcolate per ciascun ruolo neuronale (Interneuron, Motor, Muscle, Sensorimotor, Sensory).

Dall’analisi emerge una chiara organizzazione gerarchica della rete. Gli interneuroni presentano valori mediamente più elevati di Betweenness e Degree, indicando un ruolo centrale nei percorsi di comunicazione e una maggiore connettività rispetto agli altri neuroni.

La Closeness centrality risulta più alta e omogenea per interneuroni e sensorimotori, evidenziando un accesso più rapido e diretto al resto della rete. I neuroni muscle, invece, mostrano valori inferiori in quasi tutte le metriche, coerentemente con la loro posizione periferica e funzione terminale nella trasmissione del segnale.

La PageRank centrality rivela che, pur non essendo altamente connessi, i neuroni motor e muscle mantengono un’importanza funzionale elevata, poiché ricevono input da nodi centrali e influenti.

9.4 Analisi Del “Chi Parla con Chi”

nodes_df <- g_somatic %>% activate(nodes) %>% as_tibble() %>% mutate(id = row_number())
edges_df <- g_somatic %>% activate(edges) %>% as_tibble()

# Uniamo per capire Ruolo Partenza (Source) -> Ruolo Arrivo (Target)
heatmap_data <- edges_df %>%
  left_join(nodes_df %>% select(id, role_from = role), by = c("from" = "id")) %>%
  left_join(nodes_df %>% select(id, role_to = role), by = c("to" = "id")) %>%
  # Contiamo quante connessioni esistono per ogni coppia di ruoli
  group_by(role_from, role_to) %>%
  summarise(link_count = n(), .groups = 'drop')

# Dobbiamo sapere quanto sono grandi i gruppi per calcolare le connessioni possibili
group_sizes <- nodes_df %>% group_by(role) %>% summarise(N = n())

# Calcoliamo la densità
heatmap_final <- heatmap_data %>%
  left_join(group_sizes, by = c("role_from" = "role")) %>%
  rename(N_from = N) %>%
  left_join(group_sizes, by = c("role_to" = "role")) %>%
  rename(N_to = N) %>%
  mutate(
    possible_links = N_from * N_to,
    density = link_count / possible_links
  )

# visualizzazione
ggplot(heatmap_final, aes(x = role_to, y = role_from, fill = density)) +
  # Disegniamo i quadrati
  geom_tile(color = "white", size = 0.5) +
  
  # Aggiungiamo i numeri dentro (Densità in percentuale)
  geom_text(aes(label = round(density * 100, 1)), color = "black", size = 4) +
  
  # Scala Colori (Bianco -> Rosso intenso)
  scale_fill_gradient(low = "white", high = "#FF4500", name = "Densità (%)") +
  
  # Ordiniamo gli assi in modo logico (Input -> Output)
  scale_y_discrete(limits = rev(c("Sensory", "Sensorimotor", "Interneuron", "Motor", "Muscle"))) +
  scale_x_discrete(limits = c("Sensory", "Sensorimotor", "Interneuron", "Motor", "Muscle")) +
  
  labs(title = "Matrice di Densità delle Connessioni",
       subtitle = "La % indica la probabilità che un neurone X si connetta a un neurone Y.\n(Leggere: Riga -> Colonna)",
       x = "Destinazione (Target)",
       y = "Sorgente (Source)") +
  
  theme_minimal() +
  theme(axis.text = element_text(size = 11, face = "bold"),
        panel.grid = element_blank())

9.4.1 Analisi della Matrice di Densità delle Connessioni

La matrice di densità delle connessioni mostra, per ogni coppia di ruoli neuronali, la percentuale di connessioni dirette che esistono rispetto a quelle teoricamente possibili (Sorgente → Destinazione). Questo tipo di visualizzazione consente di osservare in modo sintetico la probabilità di interazione fra gruppi funzionali nella rete neurale del C. elegans.

Dall’analisi emerge che alcuni gruppi presentano una connettività interna più densa, mentre altri mostrano connessioni preferenziali verso ruoli specifici. In particolare, i neuroni Sensorimotor mostrano una densità interna relativamente alta (≈7.1 %), indicando che questi nodi tendono a connettersi frequentemente tra loro. Anche gli Interneuroni mostrano una forte connettività interna (≈6.4 %), coerente con il loro ruolo centrale di integrazione delle informazioni nella rete.

Le connessioni dai neuroni Sensory verso gli Interneuroni sono anch’esse tra le più elevate (circa 5.6 %), riflettendo il classico percorso di trasmissione dal livello sensoriale a quello integrazione. Al contrario, le connessioni dirette dai Sensory ai Motor risultano molto rare (≈1.5 %), suggerendo che il flusso di informazione percorre tipicamente uno o più strati di processi intermedi prima di raggiungere l’effettore motore.

Nel complesso, lo schema di densità delle connessioni evidenzia una struttura direzionale e gerarchica della rete: neuroni sensoriali, integrati da interneuroni e sensorimotori, e infine trasmessi ai motoneuroni che controllano i muscoli.

9.5 Analisi della connettività con l’investimento sinaptico

df_strength <- g_somatic %>%
  activate(nodes) %>%
  mutate(
    # Degree: Numero di partner unici (Quanti amici ho)
    Degree = centrality_degree(mode = "all"),
    
    # Strength: Somma totale dei pesi (Quante sinapsi totali ho)
    Strength = centrality_degree(weights = weight, mode = "all"),
    
    # Calcoliamo il rapporto: "In media, quante sinapsi dedico a ogni amico?"
    Avg_Weight_Per_Link = Strength / Degree
  ) %>%
  as_tibble() %>%
  # Rimuoviamo i nodi isolati (Degree 0) per evitare errori nei logaritmi
  filter(Degree > 0)

ggplot(df_strength, aes(x = Degree, y = Strength)) +
  # Punti: Colore = Ruolo, Dimensione = Peso Medio (per vedere gli "specialisti")
  geom_point(aes(color = role, size = Avg_Weight_Per_Link), alpha = 0.7) +
  
  # Aggiungiamo una linea che ci dice l'andamento medio
  geom_smooth(method = "lm", color = "black", linetype = "dashed", se = FALSE, size = 0.5) +
  
  # Etichette solo per i nodi speciali (quelli con connessioni fortissime)
  geom_text_repel(data = filter(df_strength, Avg_Weight_Per_Link > quantile(Avg_Weight_Per_Link, 0.95)),
                  aes(label = name), size = 3, box.padding = 0.5, max.overlaps = 20) +
  
  # Scale Logaritmiche (Fondamentali per le Power Laws)
  scale_x_log10(breaks = c(1, 2, 5, 10, 20, 50, 100)) +
  scale_y_log10(breaks = c(1, 10, 100, 1000)) +
  
  scale_color_manual(values = my_node_colors) +
  labs(title = "Relazione tra Connettività e Investimento Sinaptico",
       subtitle = "Scala Log-Log. Punti sopra la linea = Connessioni più forti della media.",
       x = "Degree (Numero di Partner)",
       y = "Strength (Numero Totale di Sinapsi)",
       size = "Peso Medio\n(Sinapsi/Link)") +
  theme_minimal() +
  theme(legend.position = "right")

## `geom_smooth()` using formula = 'y ~ x'

9.5.1 Spiegazione

Il grafico illustra la relazione tra il numero di partner sinaptici di ciascun neurone (Degree) e il numero totale di sinapsi (Strength), in scala log-log. Ogni punto rappresenta un neurone del C. elegans: il colore identifica il ruolo funzionale, mentre la dimensione del punto indica il peso medio delle connessioni (sinapsi per link). La linea tratteggiata rappresenta la relazione media di proporzionalità tra numero di connessioni e forza sinaptica totale.

Il grafico evidenzia una legge di proporzionalità diretta (relazione quasi lineare in scala log-log) tra il numero di partner (Degree) e la forza sinaptica totale (Strength).

Consolidamento del Core (Sopra la Linea): I Motoneuroni mostrano un “peso medio per link” superiore alla norma. Il sistema nervoso investe metabolicamente su di loro creando connessioni fisicamente robuste (molte sinapsi per singolo contatto) per garantire un’alta affidabilità nella trasmissione dei comandi motori.

10 Analisi Clustering

Un ulteriore domanda che mi sono posto è come la rete somatica si divide in comunità

g_undirected <- g_somatic %>%
  activate(edges) %>%
  select(from, to, weight) %>%
  as.igraph() %>% # Passiamo a igraph per la conversione sicura
  igraph::as.undirected(mode = "collapse", edge.attr.comb = list(weight = "sum")) %>%
  as_tbl_graph() 

# Calcolo delle comunità aggiungendo agli atributi

g_comm_final <- g_undirected %>%
  activate(nodes) %>%
  mutate(
    # Ora il grafo è pulitissimo, ha solo i pesi sommati
    comm_walktrap = as.factor(group_walktrap(weights = weight)),
    comm_infomap  = as.factor(group_infomap(weights = weight)), 
    comm_louvain  = as.factor(group_louvain(weights = weight))
  )

# Modularità
calc_Q_undir <- function(cluster_col) {
  # Prendiamo l'appartenenza
  memb <- g_comm_final %>% 
    activate(nodes) %>% 
    pull(!!sym(cluster_col)) %>% 
    as.numeric()
  
  # Prendiamo i pesi (dagli archi)
  weights_vector <- g_comm_final %>% 
    activate(edges) %>% 
    pull(weight)
  
  # Sicurezza
  if(any(is.na(weights_vector))) weights_vector[is.na(weights_vector)] <- 1
  
  modularity(g_comm_final, memb, weights = weights_vector)
}

results <- tibble(
  Algorithm = c("Walktrap", "Infomap", "Louvain"),
  Modularity = c(
    calc_Q_undir("comm_walktrap"),
    calc_Q_undir("comm_infomap"),
    calc_Q_undir("comm_louvain")
  ),
  Clusters = c(
    length(unique(pull(g_comm_final, comm_walktrap))),
    length(unique(pull(g_comm_final, comm_infomap))),
    length(unique(pull(g_comm_final, comm_louvain)))
  )
)

10.0.1 Bar Plot modularità

plot_bar <- ggplot(results, aes(x = reorder(Algorithm, -Modularity), y = Modularity, fill = Algorithm)) +
  geom_col(width = 0.6, alpha = 0.9, color = "black") + 
  
  # Aggiungiamo il valore numerico sopra la barra
  geom_text(aes(label = round(Modularity, 3)), vjust = -0.5, size = 5, fontface = "bold") +
  
  # Aggiungiamo il numero di cluster trovati dentro la barra (o sotto)
  geom_label(aes(y = 0.05, label = paste("Cluster:", Clusters)), 
             fill = "white", alpha = 0.8, size = 3.5) +
  
  scale_fill_brewer(palette = "Set2") +
  
  labs(title = "Qualità della Partizione Anatomica",
       subtitle = "Confronto della Modularità (Q) su rete non diretta.\nValori più alti indicano una struttura a comunità più definita.",
       y = "Modularità (Q)",
       x = NULL) +
  
  theme_minimal() +
  theme(
    legend.position = "none",
    plot.title = element_text(face = "bold", size = 16),
    axis.text.x = element_text(face = "bold", size = 12),
    panel.grid.major.x = element_blank()
  )

print(plot_bar)

10.0.2 Rappresentazione visiva dei cluster

plot_network <- ggraph(g_comm_final, layout = "fr") + # Layout Fruchterman-Reingold (fisico)
  
  geom_edge_link(alpha = 0.15, color = "grey60") +
  
  geom_mark_hull(aes(x, y, fill = comm_louvain, label = comm_louvain), 
                 concavity = 3,       # Quanto è "morbida" la forma
                 alpha = 0.2,         # Trasparenza dell'area
                 size = 0.3,          # Spessore del contorno
                 show.legend = FALSE) +
  
  geom_node_point(aes(color = comm_louvain), size = 2.5, show.legend = FALSE) +
  
  scale_fill_viridis_d(option = "plasma") + # Colori ad alto contrasto
  scale_color_viridis_d(option = "plasma") +
  
  labs(title = "Mappa dei Gangli Anatomici (Louvain)",
       subtitle = "I gruppi evidenziati rappresentano regioni ad alta densità di connessione fisica.") +
  
  theme_graph()

# Visualizza la mappa
print(plot_network)

10.1 Analisi Cluster diretto

# Usiamo g_somatic originale che contiene le frecce (direzione)
g_comm_directed <- g_somatic %>%
  activate(nodes) %>%
  mutate(
    # Infomap: Segue il flusso (Flow-based)
    comm_infomap_dir = as.factor(group_infomap(weights = weight)),
    # Walktrap: Random walk (Topology-based)
    comm_walktrap_dir = as.factor(group_walktrap(weights = weight))
  )

calc_Q_dir <- function(graph, cluster_col) {
  memb <- graph %>% activate(nodes) %>% pull(!!sym(cluster_col)) %>% as.numeric()

  weights_vector <- graph %>% activate(edges) %>% pull(weight)
  if(any(is.na(weights_vector))) weights_vector[is.na(weights_vector)] <- 1
  
  modularity(graph, memb, weights = weights_vector)
}

results_dir <- tibble(
  Algorithm = c("Infomap (Flow)", "Walktrap (Random Walk)"),
  Modularity = c(
    calc_Q_dir(g_comm_directed, "comm_infomap_dir"),
    calc_Q_dir(g_comm_directed, "comm_walktrap_dir")
  ),
  Clusters = c(
    length(unique(pull(g_comm_directed, comm_infomap_dir))),
    length(unique(pull(g_comm_directed, comm_walktrap_dir)))
  )
)

plot_bar_dir <- ggplot(results_dir, aes(x = Algorithm, y = Modularity, fill = Algorithm)) +
  geom_col(width = 0.5, alpha = 0.9, color = "black") +
  geom_text(aes(label = round(Modularity, 3)), vjust = -0.5, size = 5, fontface = "bold") +
  geom_label(aes(y = 0.05, label = paste("Cluster:", Clusters)), 
             fill = "white", alpha = 0.8, size = 3.5) +
  scale_fill_manual(values = c("#E76BF3", "#00BFC4")) + # Colori diversi dal precedente per distinguere
  labs(title = "Qualità della Partizione Funzionale (Diretta)",
       subtitle = "Confronto su come gli algoritmi catturano il flusso dell'informazione.",
       y = "Modularità Diretta (Q)", x = NULL) +
  theme_minimal() +
  theme(legend.position = "none",
        plot.title = element_text(face = "bold", size = 16))

print(plot_bar_dir)

set.seed(123)
plot_net_dir <- ggraph(g_comm_directed, layout = "fr") + 
  
  # Archi: Usiamo le frecce per ricordare che è diretto
  geom_edge_link(arrow = arrow(length = unit(2, 'mm')), 
                 end_cap = circle(2, 'mm'),
                 alpha = 0.1, color = "grey60") +
  
  # Aree dei Moduli (Infomap)
  geom_mark_hull(aes(x, y, fill = comm_infomap_dir, label = comm_infomap_dir), 
                 concavity = 3, alpha = 0.2, size = 0.3, show.legend = FALSE) +
  
  # Nodi
  geom_node_point(aes(color = comm_infomap_dir), size = 2.5, show.legend = FALSE) +
  
  scale_fill_viridis_d(option = "turbo") + # Palette "Turbo" ottima per distinguere molti gruppi
  scale_color_viridis_d(option = "turbo") +
  
  labs(title = "Mappa dei Circuiti Funzionali (Infomap)",
       subtitle = "I cluster rappresentano unità di processamento del segnale (Sorgenti -> Destinazioni).") +
  theme_graph()

print(plot_net_dir)

(Anatomia - Louvain): I neuroni sono raggruppati perché sono vicini. gruppi compatti.

(Fisiologia - Infomap): Mostra i Processi. I neuroni sono raggruppati perché si parlano. gruppi “allungati” o mescolati, perché seguono la strada del segnale (Sensore -> Interneurone).

11 Analisi Collegamenti spaziali

Un ulteriore domanda che mi sono posto è la seguente: i neuroni connessi sono fisicamente vicini? La maggior parte degli archi è corta o ci sono molti archi lunghi, e come cambia la proabilità di connetettersi all’aumentare della distanza?

nodes_df <- g_somatic %>% activate(nodes) %>% as_tibble()
edges_df <- g_somatic %>% activate(edges) %>% as_tibble()

wiring_df <- edges_df %>%
  mutate(
    # Coordinate SORGENTE (Source)
    x_s = nodes_df$x_phys[from],
    y_s = nodes_df$y_phys[from],
    z_s = nodes_df$z_phys[from],
    
    # Coordinate DESTINAZIONE (Target)
    x_t = nodes_df$x_phys[to],
    y_t = nodes_df$y_phys[to],
    z_t = nodes_df$z_phys[to]
  ) %>%
  # Calcolo Euclideo 3D
  mutate(
    Distance = sqrt((x_s - x_t)^2 + (y_s - y_t)^2 + (z_s - z_t)^2)
  ) %>%
  filter(!is.na(Distance)) # Rimuoviamo eventuali NA

# (Wiring Cost) 
p_hist <- ggplot(wiring_df, aes(x = Distance)) +
  geom_histogram(bins = 50, fill = "#2C3E50", color = "white", alpha = 0.8) +
  geom_vline(aes(xintercept = mean(Distance)), color = "red", linetype = "dashed", size = 1) +
  labs(title = "Wiring Cost: Lunghezza delle Connessioni",
       subtitle = "La maggior parte dei collegamenti è corta",
       x = "Distanza Euclidea", y = "Frequenza") +
  theme_minimal()

# Calcoliamo tutte le distanze possibili tra TUTTI i nodi (Reali + Non connessi)
coords_matrix <- nodes_df %>% 
  select(x_phys, y_phys, z_phys) %>% 
  filter(!is.na(x_phys)) %>% 
  as.matrix()

# Matrice delle distanze di tutti contro tutti
dist_matrix <- as.matrix(dist(coords_matrix))
# Prendiamo solo il triangolo superiore (valori unici)
all_possible_distances <- dist_matrix[upper.tri(dist_matrix)]

# Creiamo 30 "cestini" di distanza
bins <- seq(0, max(all_possible_distances), length.out = 30)

# Contiamo quanti archi REALI cadono in ogni cestino
h_real <- hist(wiring_df$Distance, breaks = bins, plot = FALSE)
# Contiamo quante coppie POSSIBILI cadono in ogni cestino
h_poss <- hist(all_possible_distances, breaks = bins, plot = FALSE)

prob_df <- tibble(
  Distance = h_real$mids,           
  P_Connection = h_real$counts / h_poss$counts # Probabilità
) %>%
  filter(h_poss$counts > 0) # Evitiamo divisioni per zero

p_decay <- ggplot(prob_df, aes(x = Distance, y = P_Connection)) +
  geom_point(size = 3, color = "darkred", alpha = 0.7) +
  geom_smooth(method = "loess", color = "black", linetype="dashed", se = FALSE) +
  labs(title = "Distance Decay",
       subtitle = "Probabilità di connessione vs Distanza Fisica",
       x = "Distanza", y = "Probabilità P(d)") +
  theme_minimal()

# Unione Grafici
p_hist / p_decay

## `geom_smooth()` using formula = 'y ~ x'

• Analisi del Wiring Cost: I risultati mostrano una distribuzione delle lunghezze degli archi fortemente asimmetrica a destra. La maggior parte delle connessioni sinaptiche copre distanze fisiche minime, suggerendo una forte pressione evolutiva verso la minimizzazione del volume assonale (Wiring Economy). Tuttavia, la presenza di una ‘coda lunga’ (long-tail) indica l’esistenza di costose connessioni a lungo raggio, essenziali per l’integrazione funzionale tra parti distanti distanti (es. testa-coda).

• Analisi del Distance Decay: Il grafico della probabilità di connessione rivela un chiaro effetto di ‘Distance Decay’. La probabilità che due neuroni siano connessi (\(P(d)\)) decade monotonicamente all’aumentare della loro distanza fisica. ma allo stesso tempo tende a risalire con la distanza, facendo in modo che sia pià probabile andare a trovare una connessione distante fra due neuroni piuttosto che a media distanza (Connessione distante indica un collegamento ad esempio fra testa e coda necessaria per il movimento). Andando a confrontarlo anche con la baseline globale ,ovvero prendendo due neuroni nel verme senza sapere essi dove siano, hai solo una probabilità del 1.98% che quest’ultimi siano connessi contro quasi il 10%, se sai che due neuroni sono vicini. La vicinanza fisica moltiplica la probabilità di connessione di un fattore 5x o 6x rispetto alla media casuale.

Andando a fare un esempio ,le connessioni fra neuroni sono come le amicizie umane ho tantissimi amici nel mio quartiere (costo basso), e pochissimi amici in Australia (costo alto). La probabilità che io sia amico di qualcuno crolla drasticamente man mano che ci allontaniamo da casa mia.

12 Analisi Small World

# Calcoliamo i valori veri del tuo grafo g_somatic
real_clustering <- transitivity(g_somatic, type = "global")
real_avg_path   <- mean_distance(g_somatic, directed = TRUE)

CL_ratio <- real_clustering/real_avg_path

n_nodes <- vcount(g_somatic)
n_edges <- ecount(g_somatic)
is_dir  <- is_directed(g_somatic)

set.seed(123)
sim_results <- replicate(100, {
  # Crea grafo casuale (Erdos-Renyi)
  g_rand <- sample_gnm(n = n_nodes, m = n_edges, directed = is_dir)
  
  # Calcola metriche
  c(
    C_rand = transitivity(g_rand, type = "global"),
    L_rand = mean_distance(g_rand, directed = is_dir)
  )
})

# Medie dei grafi casuali
rand_clustering <- mean(sim_results["C_rand",], na.rm = TRUE)
rand_avg_path   <- mean(sim_results["L_rand",], na.rm = TRUE)

# Gamma (Quanto siamo più clusterizzati del caso?) -> Vuoi che sia >> 1
gamma <- real_clustering / rand_clustering

# Lambda (Quanto siamo efficienti rispetto al caso?) -> Vuoi che sia ≈ 1
lambda <- real_avg_path / rand_avg_path

# Sigma (Indice Small-World) -> Se > 1, è Small-World
sigma <- gamma / lambda


# Istogramma del Clustering Locale
# Calcoliamo il clustering per ogni singolo nodo
local_clust <- transitivity(g_somatic, type = "local")
df_clust <- tibble(Clustering = local_clust) %>% filter(!is.na(Clustering))

p1 <- ggplot(df_clust, aes(x = Clustering)) +
  geom_histogram(bins = 30, fill = "#5D8AA8", color = "black", alpha = 0.8) +
  labs(title = "Distribuzione Coefficienti di Clustering",
       subtitle = paste("Media Reale:", round(real_clustering, 3), "vs Random:", round(rand_clustering, 3)),
       y = "Frequenza", x = "Clustering Coefficient") +
  theme_minimal()

#  Istogramma Lunghezza Cammini 
# Calcoliamo TUTTI i cammini minimi tra tutte le coppie (può essere lento se la rete è enorme)
dist_table <- distance_table(g_somatic)
# Ricostruiamo i dati per il plot (distance_table dà i conteggi)
paths_df <- tibble(
  Length = 1:length(dist_table$res),
  Count = as.numeric(dist_table$res)
) %>% filter(Count > 0)

p2 <- ggplot(paths_df, aes(x = Length, y = Count)) +
  geom_col(fill = "#E39E47", color = "black", alpha = 0.8, width = 1) +
  labs(title = "Distribuzione Lunghezza Cammini",
       subtitle = paste("Media Reale:", round(real_avg_path, 3), "vs Random:", round(rand_avg_path, 3)),
       y = "Numero di Percorsi", x = "Lunghezza (Salti)") +
  theme_minimal()

# Mostra i grafici
p1 / p2

• ANALISI SMALL-WORLD (Watts-Strogatz)
  • Clustering Reale: 0.2007 | Clustering Casuale: 0.0392 | Gamma: 5.13 (Target: >> 1)
  • Cammino Medio Reale: 4.8907 | Cammino Medio Casuale: 3.1856 | Lambda: 1.54 (Target: ≈ 1)
  • Rapporto C/L: 0.041
  • INDICE SIGMA: 3.3386

Il confronto che è stato eseguito con i modelli nulli (Grafi Casuali di Erdos-Renyi) ha permesso di classificare l’architettura della rete.

Metriche di Watts-Strogatz

• Clustering (\(C \gg C_{rand}\)): La rete mostra un coefficiente di clustering enormemente superiore al caso (Rapporto \(\gamma \approx 40-50\)), indicando una specializzazione locale estrema.

• Path Length (\(L \approx L_{rand}\)): La lunghezza media dei cammini è comparabile, seppur leggermente superiore, a quella casuale.

• Indice Sigma (\(\sigma \gg 1\)): Con un valore di \(\sigma > 1\), la rete è classificata teoricamente come Small-World.

Il rapporto diretto \(C/L \approx 0.041\) è inferiore ai modelli teorici ideali. Questo non indica inefficienza, ma riflette il vincolo morfologico: il fatto che’organismo impone un limite fisico alla minimizzazione del cammino medio (\(L\)) è dovuto al motivo che i neuroni sono cellule fisiche con cavi che occupano spazio e consumano energia per essere costruiti e mantenuti, che la rete compensa massimizzando l’efficienza locale (\(C\)).

In conclusione la rete preferisci collegamenti corti e locali e usa solo pochi collegamenti lunghi costosi. Questo la rende una small world accettabile, ma non perfetta

13 Analisi Max K-Core vs Rich Club e Correlazioni (Hub/Auth vs Degree)

Un ulteriore analisi che sono andato ad eseguire è stata quella di confrontare la Shell più profonda (Cricca) con i neuroni che hanno tanti contatti, per vedere se c’è sovrapposizione.

g_check <- g_somatic %>%
  activate(nodes) %>%
  mutate(
    degree_total = centrality_degree(mode = "all"),
    k_core_index = node_coreness(mode = "all")
  ) %>%
  as_tibble()
max_k_val <- max(g_check$k_core_index)

top15_list <- g_check %>%
  arrange(desc(degree_total)) %>%
  slice(1:15) %>%
  select(name, degree_total, k_core_index) %>%
  mutate(
    Status = ifelse(k_core_index == max_k_val, "CORE HUB", "bridge CONNECTOR"),
    Interpretation = ifelse(k_core_index == max_k_val, 
                            "Membro dell'Elite (Rich + Core)", 
                            "Ponte verso la Periferia")
  )

# Tabella Dettagliata
top15_list %>%
  kable(format = "markdown", 
        col.names = c("Neurone", "Grado Totale", "K-Shell Index", "Ruolo", "Interpretazione"),
        align = "lcccl",
        caption = "Analisi Topologica dei Top 15 Hub: Core vs Connector")

Analisi Topologica dei Top 15 Hub: Core vs Connector

Neurone	Grado Totale	K-Shell Index	Ruolo	Interpretazione
AVAR	98	12	CORE HUB	Membro dell’Elite (Rich + Core)
AVAL	90	12	CORE HUB	Membro dell’Elite (Rich + Core)
AVBL	60	12	CORE HUB	Membro dell’Elite (Rich + Core)
PVCL	59	12	CORE HUB	Membro dell’Elite (Rich + Core)
PVCR	58	12	CORE HUB	Membro dell’Elite (Rich + Core)
AVDR	57	12	CORE HUB	Membro dell’Elite (Rich + Core)
DVA	54	11	bridge CONNECTOR	Ponte verso la Periferia
AVBR	53	12	CORE HUB	Membro dell’Elite (Rich + Core)
AVEL	52	12	CORE HUB	Membro dell’Elite (Rich + Core)
AVER	51	11	bridge CONNECTOR	Ponte verso la Periferia
AVDL	46	12	CORE HUB	Membro dell’Elite (Rich + Core)
RIAR	44	11	bridge CONNECTOR	Ponte verso la Periferia
RIAL	42	11	bridge CONNECTOR	Ponte verso la Periferia
HSNR	41	12	CORE HUB	Membro dell’Elite (Rich + Core)
RIMR	38	11	bridge CONNECTOR	Ponte verso la Periferia

# Riassunto Finale
n_core <- sum(top15_list$k_core_index == max_k_val)
n_conn <- 15 - n_core

cat("Su 15 Hub analizzati:\n")

## Su 15 Hub analizzati:

cat("-", n_core, "sono CORE HUBS (Vivono nel centro profondo)\n")

## - 10 sono CORE HUBS (Vivono nel centro profondo)

cat("-", n_conn, "sono CONNECTOR HUBS (Vivono in gusci più esterni)\n")

## - 5 sono CONNECTOR HUBS (Vivono in gusci più esterni)

Analizzando la tabella possiamo notare che su 21 neuroni appartenenti al Max K-shell (12) 10 di questi appartengono al rich club, sui 15 totali.

La maggioranza degli hubs risiede nel max K-Core, quasi il 66%.

Ma l’analisi più interessante è che la sovrapposizione non è completa. 5 dei 15 neuroni appartenti al rich club non fanno parte della Max K-shell. Questi nodi si possono interpetare come dei ponti verso l’esterno del sistema (come ad esempio trasferire un input al muscolo).

Inoltre 11 dei neuroni presenti all’interno del max core non sono top hubs, questo implica che pur avendo un numero moderato di connessioni, garantiscono la resilienza del nucleo pure se gli hub venissero rimossi.

13.1 Grafico rappresentativo del Max Core e Overlap Rich-Club

g_core_viz <- g_somatic %>%
  activate(nodes) %>%
  mutate(
    # Calcoliamo la coreness (indice del guscio)
    k_shell = node_coreness(mode = "all"),
    
    # Identifichiamo il valore massimo
    max_k_val = max(k_shell),
    
    # Flag: Sei nel nucleo massimo?
    is_max_core = k_shell == max_k_val,
    
    # Etichetta
    core_label = ifelse(is_max_core, name, NA),
    
    node_col = ifelse(is_max_core, "red", "grey90"),
    node_alpha = ifelse(is_max_core, 1, 0.2),
    node_size = ifelse(is_max_core, 5, 2)
  ) %>%
  # Passiamo agli archi
  activate(edges) %>%
  mutate(
    # Qui usiamo .N() per guardare le proprietà dei nodi mentre siamo negli archi
    from_core = .N()$is_max_core[from],
    to_core = .N()$is_max_core[to],
    
    # Un arco è "Core" solo se connette due nodi del nucleo
    is_core_edge = from_core & to_core,
    edge_col = ifelse(is_core_edge, "red", "grey95"),
    edge_alpha = ifelse(is_core_edge, 0.6, 0.05)
  )

max_k_shell_value <- g_core_viz %>% 
  activate(nodes) %>% 
  pull(k_shell) %>% 
  max()

set.seed(123) 
ggraph(g_core_viz, layout = "fr") + 
  
  geom_edge_link(aes(filter = !is_core_edge), color = "grey95", alpha = 0.1) +
  
  geom_edge_link(aes(filter = is_core_edge), color = "#E74C3C", alpha = 0.6, width = 0.8) +
  
  geom_node_point(aes(filter = !is_max_core), color = "grey80", size = 2, alpha = 0.2) +
  
  geom_node_point(aes(filter = is_max_core), color = "#E74C3C", size = 5) +
  
  geom_node_text(aes(label = core_label), repel = TRUE, 
                 size = 3, fontface = "bold", color = "black",
                 max.overlaps = 20) +
  
  labs(title = paste("Visualizzazione del MAX K-CORE (K =", max_k_shell_value, ")"),
       subtitle = "I nodi rossi rappresentano il nucleo strutturale più profondo della rete.\nIl resto del sistema nervoso è in grigio.",
       caption = "Nota: Questo grafico isola la 'Cipolla' centrale della rete (Shell massima).") +
  
  theme_graph()

g_stats <- g_somatic %>%
  activate(nodes) %>%
  mutate(
    in_degree = centrality_degree(mode = "in"),
    out_degree = centrality_degree(mode = "out"),
    total_degree = centrality_degree(mode = "total"),
    pagerank = centrality_pagerank(),
    k_core = node_coreness(mode = "all"),
    hub_score = centrality_hub(),
    auth_score = centrality_authority(),
    betweenness = centrality_betweenness(directed = TRUE) 
  ) %>%
  as_tibble()

# Definiamo le soglie
degree_threshold <- quantile(g_stats$total_degree, 0.90) # Top 10% Grado
max_core_val <- max(g_stats$k_core)

g_overlap <- g_stats %>%
  mutate(
    Is_Rich = total_degree >= degree_threshold,
    Is_MaxCore = k_core == max_core_val,
    Category = case_when(
      Is_Rich & Is_MaxCore ~ "Rich Club + Max Core (L'Elite)",
      Is_Rich & !Is_MaxCore ~ "Solo Rich Club (Hub Periferici)",
      !Is_Rich & Is_MaxCore ~ "Solo Max Core (Deep Connectors)",
      TRUE ~ "Altri"
    )
  )

top_betweenness_names <- g_overlap %>%
  arrange(desc(betweenness)) %>% # Ordina dal più alto al più basso
  slice(1:10) %>%                # Prendi i primi 10
  pull(name)                # Estrai solo i nomi

# Applichiamo l'etichetta solo se il nome è nella lista dei Top 10
g_overlap <- g_overlap %>%
  mutate(label_to_show = ifelse(name %in% top_betweenness_names, name, NA))

p_overlap <- ggplot(g_overlap, aes(x = total_degree, y = k_core)) +

  geom_jitter(aes(color = Category, size = betweenness), 
              alpha = 0.7, width = 0.2, height = 0.2)+
  
  geom_text_repel(aes(label = label_to_show),
                  box.padding = 0.6,
                  point.padding = 0.3,
                  max.overlaps = Inf,
                  size = 4,             # Un po' più grande per leggibilità
                  fontface = "bold",
                  color = "black") +
  
  geom_vline(xintercept = degree_threshold, linetype = "dashed", color = "grey50") +
  geom_hline(yintercept = max_core_val - 0.5, linetype = "dashed", color = "grey50") +
  
  scale_color_manual(values = c(
    "Altri" = "grey85",
    "Rich Club + Max Core (L'Elite)" = "#E74C3C",  # Rosso
    "Solo Rich Club (Hub Periferici)" = "#F39C12", # Arancione
    "Solo Max Core (Deep Connectors)" = "#3498DB"  # Blu
  )) +
  
  # Gestione Dimensioni (Range controlla quanto grandi diventano le bolle)
  scale_size_continuous(range = c(2, 10), name = "Betweenness Centrality") +
  
  labs(title = "Structural Overlap: Rich Club vs K-Core",
       subtitle = "Confronto tra i nodi più connessi (Rich) e i nodi più centrali (Core)",
       caption = "I punti rossi rappresentano i neuroni che sono SIA Hub che parte del Nucleo.",
       x = "Grado Totale (k)", y = "Indice K-Shell (Coreness)") +
  theme_minimal() +
  theme(legend.position = "bottom",
        legend.box = "vertical")

print(p_overlap)

Il grafico scatter plot mette in relazione la quantità di connessioni (Grado Totale, asse X) con la profondità strutturale (Indice K-Shell, asse Y), calcolata considerando la connettività bidirezionale (mode=“all”).

Dal grafico emergono tre distinte popolazioni funzionali:

• L’Élite (Rich Club + Max Core): I nodi in rosso rappresentano il centro di comando della rete. Essi possiedono simultaneamente un altissimo grado (\(k > 28\)) e la massima Coreness (\(K_{max} = 12\)).

• Hub Periferici (Connector Hubs) : I nodi in arancione sono neuroni ad alto grado che tuttavia non appartengono al nucleo più profondo. Pur avendo molte connessioni, la loro posizione topologica è leggermente più esterna (\(K \approx 10-11\)). Il loro ruolo è quello di ponti informativi: raccolgono segnali diffusi dalla periferia sensoriale e li convogliano verso il centro.

• Deep Connectors (La “Colla” del Nucleo) : I nodi in azzurro (sinistra, sopra la linea tratteggiata) sono la componente strutturale più interessante. Pur non essendo “famosi” in termini di grado (hanno meno di 28 connessioni), risiedono nel guscio massimo (\(K=12\)). Questi neuroni agiscono come “connettori profondi”, garantendo la coesione interna del nucleo e impedendone la frammentazione qualora gli Hub principali venissero rimossi.

14 Correlazioni e Analisi Hub e Authority

Un ulteriore analisi che sono andato ad effettuare sono le Correlazioni (Hub/Auth vs Degree).

Inoltre dalle analisi relativi alle metriche di hub score e authority score si possono rilevare molteplici conclusioni. L’Hub Score e l’Authority Score sono due metriche di centralità derivate dall’algoritmo HITS.

• Authority Score (L’Autorevolezza - Focus sull’Input): Un nodo ha un alto punteggio di Authority se riceve connessioni da molti Hub importanti.

• Hub Score (Il Collettore - Focus sull’Output):Un nodo ha un alto punteggio di Hub se invia connessioni verso molte Authority importanti. Quindi è presente una relazione ricosiva: un buon Hub è tale perché punta a buone Authority; una buona Authority è tale perché è puntata da buoni Hub.

Inoltre il grafico Authority vs Hub Score, rappresenta esattemente un sistema biologico. i nodi in alto a destra hanno punteggi altri su entrambi gli assi, ovvero che ricevono informazioni (Authority), le elaborano e le ritrasmettono (Hub). Non sono specializzati, gestiscono l’intero flusso di informazioni Altri nodi sono sbilanciati verso uno dei due assi, indicando una specializzazione (solo ricevitori o solo trasmettitori).

il prestigio (PageRank) è dettato dalla popolarità (In-Degree), con una correlazione molto alta (0.901). questi nodi rappresentano molto probabilmente i terminali decisionali del sistema

15 Analisi di Resilienza con diverse strategie

# attacchi STATICI 
perform_attack_static <- function(graph, node_order, attack_name) {
  n_total <- vcount(graph)
  results <- numeric(length(node_order))
  g_temp <- graph
  
  for (i in seq_along(node_order)) {
    g_temp <- delete_vertices(g_temp, node_order[i])
    if (vcount(g_temp) == 0) results[i] <- 0
    else results[i] <- max(components(g_temp)$csize)
  }
  
  tibble(
    Fraction_Removed = seq_along(node_order) / n_total,
    LCC_Size = results,
    Strategy = attack_name
  )
}

# attacchi DINAMICI 
# Questa funzione accetta una metrica ("degree", "betweenness", "pagerank") e ricalcola dinamicamente ogni metrica quando si toglie un nodo
perform_attack_dynamic <- function(graph, metric_type, attack_name) {
  g_temp <- graph
  n_total <- vcount(g_temp)
  results <- numeric(n_total)
  
  
  for (i in 1:n_total) {
    if (vcount(g_temp) == 0) {
      results[i:n_total] <- 0
      break
    }
    
    # Ricalcolo della metrica sul grafo sopravvissuto
    if (metric_type == "degree") {
      vals <- degree(g_temp)
    } else if (metric_type == "betweenness") {
      vals <- betweenness(g_temp) # Lento ma letale
    } else if (metric_type == "pagerank") {
      vals <- page_rank(g_temp)$vector
    }
    
    # Trova il nodo con valore massimo (gestione spareggi: prende il primo)
    target <- names(vals)[which.max(vals)]
    
    # Rimuovi
    g_temp <- delete_vertices(g_temp, target)
    
    # Misura
    if (vcount(g_temp) == 0) results[i] <- 0
    else results[i] <- max(components(g_temp)$csize)
  }
  
  tibble(
    Fraction_Removed = (1:n_total) / n_total,
    LCC_Size = results,
    Strategy = attack_name
  )
}

# Assumiamo g_somatic caricato
g_base <- as.igraph(g_somatic)

# Liste Statiche Standard ---
set.seed(123)
target_random <- sample(V(g_base)$name)
target_degree <- V(g_base)$name[order(degree(g_base), decreasing = TRUE)]
target_betw   <- V(g_base)$name[order(betweenness(g_base), decreasing = TRUE)]

# Attacco Biologico "Smart
global_betw <- betweenness(g_base)
inter_names <- V(g_base)$name[V(g_base)$role == "Interneuron"] # Verifica se 'role' esiste
inter_betw <- global_betw[inter_names]
inter_sorted <- names(sort(inter_betw, decreasing = TRUE))
others <- setdiff(V(g_base)$name, inter_sorted)
target_inter_smart <- c(inter_sorted, sample(others))

# Calcoliamo la "Coreness" (A quale guscio appartiene ogni nodo?)
k_shells <- coreness(g_base)

# Creiamo un dataframe per ordinare i bersagli
df_kcore <- tibble(
  name = V(g_base)$name,
  k_shell = k_shells,           # Il guscio (es. 1, 2... 10)
  betw = global_betw            # La Betweenness globale
)

# Logica di Ordinamento:
# PRIMA: Ordina per K-Shell decrescente (Partiamo dal nucleo più denso, es. K=10)
# POI: A parità di guscio, uccidi chi ha più Betweenness
target_kcore_smart <- df_kcore %>%
  arrange(desc(k_shell), desc(betw)) %>%
  pull(name)

# Eseguiamo gli statici
df_static <- bind_rows(
  perform_attack_static(g_base, target_random,      "Random (Failure)"),
  perform_attack_static(g_base, target_degree,      "Degree (Static)"),
  perform_attack_static(g_base, target_inter_smart, "Interneurons (Smart Bio)"),
  perform_attack_static(g_base, target_kcore_smart, "K-Core Smart (Core+Betw)")
)

# Eseguiamo i dinamici 
df_dyn_deg  <- perform_attack_dynamic(g_base, "degree",      "Dynamic Degree")
df_dyn_betw <- perform_attack_dynamic(g_base, "betweenness", "Dynamic Betweenness") 
df_dyn_page <- perform_attack_dynamic(g_base, "pagerank",    "Dynamic PageRank")

# Uniamo tutto
data_final <- bind_rows(df_static, df_dyn_deg, df_dyn_betw, df_dyn_page)

ggplot(data_final, aes(x = Fraction_Removed, y = LCC_Size, color = Strategy)) +
  geom_line(size = 1.0, alpha = 0.8) +
  
  # Linea di riferimento 50%
  geom_hline(yintercept = vcount(g_base)/2, linetype = "dashed", color = "grey") +
  
  scale_color_manual(values = c(
    "Random (Failure)" = "#2ECC71",         # Verde (Base)
    "Interneurons (Smart Bio)" = "#FF00FF", # Magenta (Bio)
    
    "Degree (Static)" = "red",          # Rosso Chiaro
    
    "K-Core Smart (Core+Betw)" = "yellow",   # giallo 
    
    "Dynamic Degree" = "grey30",            # Grigio Scuro
    "Dynamic PageRank" = "#3498DB",         # Blu (Ottimo per diretti)
    "Dynamic Betweenness" = "black"         # Nero (La morte nera)
  )) +
  
  labs(title = "Analisi Avanzata di Resilienza",
       subtitle = "Confronto tra Statici, Ibridi e Dinamici",
       x = "Frazione di Nodi Rimossi (f)",
       y = "Dimensione Componente Gigante") +
  
  theme_minimal() +
  theme(legend.position = "right")

Come il grafico evidenza la rete è resiliente riguardo i guasti, lo si può dedurre dal fatto che la strategia casuale è la peggiore delle strategie presenti. Dall’altro canto la rete somatica è molto sensibile agli attacchi mirati, questa caratteristica è dovuta alla presenza di super hub, infatti la rimozione per grado e PageRank sono risultate le strategie più efficaci.

Questo conferma il fatto che la rete è scale free, il fatto che attaccando la rete a caso, il sistema non collassa. Questo indica la robustezza della rete. Al contrario la fragilità della rete è data dagli attacchi mirati, se si attacca una piccola frazione di nodi specificie, e li elimini, il sistema muore.

é importante sottolineare come l’attaco basato sulla degree statica (non la ricalcoliamo sempre) è già sorpendetemente efficace. Gli hub sono già noti dall’inizio

16 Conclusioni

L’analisi della rete somatica del verme C ELEGANS ha messo in luce moltepici aspeti significativi caratteristici di una rete complessa

• La rete è scale free vincolata fisicamente, supportando il fallimento dei nodi, mostrando la resilienza dei nodi.
• La rete è orientata alla creazione di Hub (Power Law), ma è vincolata alla creazione di quest’ultimi per limiit fisici (cut-off esponenziale)
• La rete mostra significative proprietà Small-World, quindi con Cammino Medio breve(soggetto a vincoli fisici) + Clustering alto.
• La formazione di comunità rispetto alla fisica reale del verme
• La probabilità che due neuroni siano connessi (\(P(d)\)) decade monotonicamente all’aumentare della loro distanza fisica.